### Make sure all the packages are installed
if (!requireNamespace("Seurat", quietly = TRUE))
install.packages("Seurat")
if (!requireNamespace("tidyverse", quietly = TRUE))
install.packages("tidyverse")
if (!requireNamespace("colorBlindness", quietly = TRUE))
install.packages("colorBlindness")
if (!requireNamespace("RColorBrewer", quietly = TRUE))
install.packages("RColorBrewer")
if (!requireNamespace("DT", quietly = TRUE))
install.packages("DT")
if (!requireNamespace("presto", quietly = TRUE))
devtools::install_github("immunogenomics/presto")
### Load all the necessary libraries
library(Seurat)
library(tidyverse)
library(colorBlindness)
library(RColorBrewer)
library(DT)7 - Differential Gene Expression & Level 1 annotation
Introduction
In this notebook we are going to look at how to interpret and visualize gene-level statistics obtained from differential expression analysis. We are not going to go into which method should be used to carry out differential gene expression analysis but we highly recommend giving a read to A comparison of marker gene selection methods for single-cell RNA sequencing data by Jeffrey M. Pullin & Davis J. McCarthy if you’re interested in digging deeper!
Some other interesting papers and twitter discussions can be found here:
Why Seurat and Scanpy’s log fold change calculations are discordant - https://twitter.com/lpachter/status/1694387749967847874?s=46.
Discrepancies between Seurat and Scanpy’s logFC - https://twitter.com/slavov_n/status/1582347828818456576
Differences in wilcoxon rank sum test p-value calculations between Seurat and Scanpy - https://twitter.com/Sanbomics/status/1693995213298266515
A comparison of marker gene selection methods for single-cell RNA sequencing data - “Overall, our results suggest that methods based on logistic regression, Student’s t-test and the Wilcoxon rank-sum test all have strong performance.”
Do you really understand log2Fold change in single-cell RNAseq data?
Key Takeaways
To annotate our clusters we need to determine which genes are differentially expressed in each one.
We can quantify these differentially expressed genes using effect size and discriminatory power metrics such as log2FC and AUC.
Differential gene expression metrics vary depending on the groups of cells we are comparing.
P values obtained from carrying out DGE analysis between clusters are inflated and should not be used.
Libraries
Load data
We’re going to be working with a dataset from the paper - Immunophenotyping of COVID-19 and influenza highlights the role of type I interferons in development of severe COVID-19 Download data from cellxgene portal.
# Download the data in data/ directory
# download.file(
# url = "https://datasets.cellxgene.cziscience.com/d8e35450-de43-451a-9979-276eac688bce.rds",
# destfile = "../data/workshop-data.rds",
# method = "wget",
# extra = "-r -p --random-wait")
# We can also use the CLI with the wget command below
# wget https://datasets.cellxgene.cziscience.com/d8e35450-de43-451a-9979-276eac688bce.rds
se <- readRDS("../data/d8e35450-de43-451a-9979-276eac688bce.rds")Generate a color palette for our cell types
# https://www.datanovia.com/en/blog/easy-way-to-expand-color-palettes-in-r/
# nb.cols <- length(unique(se$Celltype))
# mycolors <- colorRampPalette(paletteMartin)(nb.cols)
pal <- paletteMartin
names(pal) <- sort(unique(se$Celltype))Analysis
Convert ENSEMBL IDs to Gene Symbols
Right away we can see how ensembl ids are used in the rownames. Let’s transform them into their matched symbols to make them human-readable:
head(rownames(se))[1] "ENSG00000000003" "ENSG00000000005" "ENSG00000000419" "ENSG00000000457" "ENSG00000000460" "ENSG00000000938"Convert to gene symbols
gene_df <- readr::read_csv(file = "../data/cov_flu_gene_names_table.csv")
symbol_id <- data.frame(index = rownames(se)) %>%
left_join(gene_df, by = "index") %>%
pull(feature_name)
# re-create seurat object
mtx <- se@assays$RNA$data
rownames(mtx) <- symbol_id
se <- CreateSeuratObject(counts = mtx, meta.data = se@meta.data)
rm(mtx); gc() used (Mb) gc trigger (Mb) limit (Mb) max used (Mb)
Ncells 3715605 198.5 5550351 296.5 NA 5550351 296.5
Vcells 180966095 1380.7 781825193 5964.9 98304 818857589 6247.4Quick processing
se <- se %>%
NormalizeData(verbose = FALSE) %>%
FindVariableFeatures(
method = "vst",
nfeatures = 3000,
verbose = FALSE) %>%
ScaleData(verbose = FALSE, features = VariableFeatures(.)) %>%
RunPCA(verbose = FALSE)
ElbowPlot(se, ndims = 50)
Let’s run FindNeighbors and FindClusters to label our data:
se <- FindNeighbors(se, reduction = "pca") %>%
FindClusters(resolution = c(0.01, 0.05, 0.1, 0.25))Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 59572
Number of edges: 1894745
Running Louvain algorithm...
Maximum modularity in 10 random starts: 0.9956
Number of communities: 4
Elapsed time: 9 seconds
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 59572
Number of edges: 1894745
Running Louvain algorithm...
Maximum modularity in 10 random starts: 0.9829
Number of communities: 7
Elapsed time: 10 seconds
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 59572
Number of edges: 1894745
Running Louvain algorithm...
Maximum modularity in 10 random starts: 0.9721
Number of communities: 9
Elapsed time: 10 seconds
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 59572
Number of edges: 1894745
Running Louvain algorithm...
Maximum modularity in 10 random starts: 0.9532
Number of communities: 15
Elapsed time: 11 secondsLet’s compute the UMAP
se <- se %>%
RunUMAP(dims = 1:30, verbose = FALSE)We can now take a look on the UMAP space how the clusters look
se$sample_id <- se$`Sample ID`
DimPlot(
se,
group.by = c(
"RNA_snn_res.0.01", "RNA_snn_res.0.05",
"RNA_snn_res.0.1", "RNA_snn_res.0.25"),
label = TRUE)
dim_plt <- DimPlot(
se,
group.by = c("RNA_snn_res.0.05"),
label = TRUE)And the original cell type labels + the sample IDs
DimPlot(
se,
group.by = c("Celltype", "sample_id"),
label = FALSE)
For the purpose of this tutorial we’re going to go forward with resolution 0.05!
DGE Wilcoxon
The different implementations Seurat incorporates provides in FindAllMarkers compare the gene expression between 2 groups of cells. This one vs all strategy is very quick and returns the avg_log2FC. This avg_log2FC is computed as detailed here & here. Since we’re working with normalized data the log2FC can be directly computed by subtracting the average expression between both groups - log(\frac{exp1}{exp2})=log(exp1)-log(exp2)
Idents(se) <- "RNA_snn_res.0.05"
mgs <- FindAllMarkers(
se,
test.use = "wilcox",
slot = "data",
only.pos = TRUE,
logfc.threshold = 0.5,
min.pct = 0.25)Look at the results in a dynamic table:
DT::datatable(mgs, filter = "top")Look at the results in a heatmap
top10 <- mgs %>%
group_by(cluster) %>%
dplyr::filter(avg_log2FC > 1) %>%
slice_head(n = 10) %>%
ungroup()
DoHeatmap(se, features = top10$gene) + NoLegend()
Annotation
Cluster 0 & 4
Let’s look at genes that are differentially expressed
FeaturePlot(
se,
features = c("CD3D", "CD3D", "TRAC", "TRBC2", "CD8B", "CD4")) +
dim_plt
VlnPlot(
se,
features = c("CD3D", "CD3D", "TRAC", "TRBC2", "CD8B", "CD4"),
group.by = "RNA_snn_res.0.05") +
dim_plt
Clusters 0 & 4 seem to have a lot of expression of T cell related genes so at this level 1 we are going to label them as T cells.
Cluster 1
Let’s look at genes that are differentially expressed
FeaturePlot(
se,
features = c("CD14", "S100A8", "VCAN", "LYZ", "MS4A6A")) +
dim_plt
VlnPlot(
se,
features = c("CD14", "S100A8", "VCAN", "LYZ", "MS4A6A"),
group.by = "RNA_snn_res.0.05") +
dim_plt
Cluster 1 is expressing a lot of monocyte-like genes, at this level 1 we are going to label them as monocytes.
Cluster 2
Let’s look at genes that are differentially expressed
FeaturePlot(
se,
features = c("MS4A1", "CD79A", "CD79B", "IGHD", "IGHM")) +
dim_plt